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Studying spin-glass physics through analyzing their ground-state properties has a long history. 
Although there exist polynomial-time algorithms for the two-dimensional planar case, where the 
problem of finding ground states is transformed to a minimum-weight perfect matching problem, 
the reachable system sizes have been limited both by the needed CPU time and by memory re- 
quirements. In this work, we present an algorithm for the calculation of exact ground states for 
two-dimensional Ising spin glasses with free boundary conditions in at least one direction. The algo- 
rithmic foundations of the method date back to the work of Kasteleyn from the 1960s for computing 
the complete partition function of the Ising model. Using Kasteleyn cities, we calculate exact ground 
states for huge two-dimensional planar Ising spin-glass lattices (up to 3000 2 spins) within reasonable 
time. According to our knowledge, these are the largest sizes currently available. Kasteleyn cities 
were recently also used by Thomas and Middleton in the context of extended ground states on the 
torus. Moreover, they show that the method can also be used for computing ground states of planar 
graphs. Furthermore, we point out that the correctness of heuristically computed ground states can 
easily be verified. Finally, we evaluate the solution quality of heuristic variants of the Bieche et al. 
approach. 

Keywords: Ising spin glass, Kasteleyn cities, perfect matching 
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I. INTRODUCTION 

The determination of spin-glass ground states has 
raised the interest of both physicists and computer scien- 
tists. For an introduction we refer to [1-3]. On the one 
hand, an analysis of the ground-state properties sheds 
light on the ruling physics of the system. On the other 
hand, several different algorithms have been developed 
and used for the ground-state determination of different 
models. 

For the two-dimensional Edwards- Anderson model [4] 
(EA) with free boundaries in at least one direction, 
ground states can be determined exactly with fast algo- 
rithms. In fact, the problem is solvable in time bounded 
by a polynomial in the size of the input. The latter can be 
achieved by a transformation to a well known graph the- 
oretical problem — the minimum-weight matching prob- 
lem, for which efficient implementations exist. For gen- 
eral non-planar or three- or higher-dimensional lattices, 
however, calculating exact ground states is NP-hard [5]. 
Loosely speaking, this means we cannot expect to be able 
to design a polynomial-time solution algorithm. In prac- 
tice, one can use, e.g. branch-and-cut algorithms [6]. 

In this work, we focus on the polynomially solvable 
case of two-dimensional lattices with free boundaries in 
at least one direction. We first review and compare the 
main known approaches which are those of Bieche et al. 
[7] and of Barahona [8, 9]. Then we present the approach 
inspired by Kasteleyn [10]. All method basically follow 
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the same idea: An associated graph is constructed in 
which a minimum-weight perfect matching is determined 
that is used to construct an exact ground state. Differ- 
ences occur in the constructed associated graph. It turns 
out that the approach inspired by Kasteleyn is the most 
favorable. In fact, using the latter method, we can de- 
termine exact ground states for lattice sizes up to 3000 2 , 
while the possible sizes computed earlier with heuristic 
variants of the approach of Bieche et al. were consider- 
ably smaller. In a forthcoming article [11], we will an- 
alyze the physics of the system. Kasteleyn cities were 
recently also used by Thomas and Middleton [12]. While 
focussing on extended ground states on the torus, they 
show that the Kasteleyn-city approach can be success- 
fully used in the planar case, too. Furthermore, they 
compared their implementation with an implementation 
of Barahona's method. It turned out that the approach 
with Kasteleyn cities is less memory consumptive and 
faster. Apart from this recent work, we are not aware of 
other computational studies using Barahona's method. 



We show how to either prove correctness of heuristi- 
cally determined ground states or how to correct them 
using linear programming. Despite the fact that this is 
fast, it is still advantageous to use the method based on 
Kasteleyn cities. Finally, we evaluate the quality of the 
solutions generated with heuristic variants of the Bieche 
et al. approach. 



The outline of the article is as follows. In Section II, 
we introduce the model. In Section III we introduce def- 
initions necessary for the literature review in Section IV. 
Finally, we report the results in Section V. 



2 



II. THE MODEL 

In the Edwards-Anderson model N spins are placed 
on a lattice. We focus on quadratic (N = L 2 ) lattices 
with free boundary conditions in at least one direction. 
Toric boundary conditions may be applied to at most one 
lattice axis. The Hamiltonian of the system is 

^ = - JijSiSj, (i) 
<ij> 

where the sum runs over all nearest-neighbor sites. Each 
spin Si is a dynamical variable which has two allowed 
states, +1 and —I. The coupling strengths Jij between 
spins i and j are independent identically distributed ran- 
dom variables following some probability distribution. 
The concentration of anti-ferromagnetic (Jij < 0) and 
ferromagnetic bonds (Jy > 0) depends on the underly- 
ing distribution. The Gaussian and the bimodal ± J dis- 
tributions are often used. A spin configuration attaining 
the global minimum of the energy function TL is called a 
ground state. 

III. PRELIMINARIES 

In this section, we briefly summarize some basic defini- 
tions from graph theory. For further details, we refer to 
[13-15] and the references therein. We associate a spin- 
glass instance with a graph G — (V, E) with vertices V 
(spin sites) and edges E (bonds). The edge set consists of 
unordered pairs (i, j), with i, j G V. More specifically, for 
two-dimensional K x L lattices with free boundaries, the 
graph is called a grid graph and is denoted by Gk.l- In 
case periodic boundaries are present in one direction, we 
call the graph a half-torus or a ring. The degree deg(v) 
of a vertex v is the number of edges (v,Wi) G E inci- 
dent at v. A path, n = v\, V2, Vk, Vi G V, is a se- 
quence of vertices such that (vi, v^), (v2, V3), (vfc-i , 
are edges of G and the Vi are distinct. A closed path 
n = Vi,V2, Vk, Vi is a cycle. 

In many applications a rational cost or a weight w(e) is 
associated with an edge e. Let G — (V, E) be a weighted 
graph. For each (possibly empty) subset Q C V, a cut 
6(Q) in G is the set of all edges with one vertex in Q 
and the other in V \ Q. The weight of a cut is given 
by w(S(Qj) = w(e). A minimum cut (min- 

cut) asks for a cut 5(Q) with minimum weight among 
all vertex sets Q C V. Let K n denote the complete graph 
with n vertices and edge set E = V x V. A subgraph 
of G is a graph Gh such that every vertex of Gh is a 
vertex of G, and every edge of Gh is an edge in G also. 
G = (V, E) is called Eulerian if and only if each vertex of 
G has even degree. A graph G is planar if it can be drawn 
in the plane in such a way that no two edges meet each 
other except at a vertex. Any such drawing is called a 
planar drawing. Any planar drawing of a graph G divides 
the plane into regions, called faces. One of these faces is 



unbounded, and called the outer face or unbounded face. 
A geometric dual graph [16] Gd of a connected planar 
graph G is a graph Gd with the property that it has 
a vertex for each face of G and an edge for each edge 
touching two neighboring faces in G. 

A matching in a graph G = (V, E) is a set of edges 
M C E such that no vertex of G is incident with more 
than one edge in M . A matching M is perfect if ev- 
ery vertex is incident to an edge in the matching. A 
maximum matching is a matching of maximum weight 
w(M) = w(e). Solving the perfect matching prob- 

eGM 

lem on general graphs in time bounded by a polynomial in 
the size of the input remained an elusive goal for a long 
time until Edmonds [17, 18] gave the first polynomial- 
time algorithm — the blossom algorithm. More details 
about matching theory can be found in [19]. 



IV. REVIEW OF THE KNOWN ALGORITHMIC 
APPROACHES 

Bicche et al. [7] showed that the problem of finding a 
ground state for two-dimensional planar Ising spin glasses 
can be transformed to a well known graph theoretical 
problem — the minimum-weight perfect matching prob- 
lem (MWPM) on general graphs. The method follows 
the scheme shown in Algorithm 1 in which an optimum 
matching is used to construct a spin configuration mini- 
mizing the total energy. 

Most commonly used exact methods, like the ap- 
proaches of Bieche et al. [7] and Barahona [8, 9], follow 
this scheme. In the following, we briefly summarize these 
two methods. Afterwards, we present a method follow- 
ing the construction introduced by Kasteleyn [10]. More 
details can be found in the recent tutorial [1] on algo- 
rithms for computing ground states in two-dimensional 
Ising spin glasses. 



Review of Exact Methods 

1. The Approach of Bieche et al. 

Bieche et al. [7] consider the weighted grid graph 
Gk,l = (V, E) where each vertex i G V is assigned an 
initial spin value Sf = ±1. Each edge e = receives 
a weight w(e) = —JijSfSj, cf. Fig. 1. Often, the trivial 
configuration S° = +1 VSi, i G V is used. 

An instance can not only be described in terms of spins 
and bonds, but also by frustrated plaquettes and paths 
of broken edges. Plaquettes consist of the 4-cycles in the 
graph. An edge is said to be satisfied if it attains its min- 
imal weight (—JijSfSj = — I |), otherwise it is called 
unsatisfied. A plaquette is frustrated if there is no spin 
configuration satisfying all edges. In this case the pla- 
quette has an odd number of negative edges. For the 
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Algorithm 1 Calculate a ground state of a K x L spin glass 
Input: Planar grid graph Gk,l 

Output: Spin configuration S minimizing the total energy H 

1. Construct an appropriate dual graph G 

2. Calculate a minimum-weight perfect matching M in G 

3. Use M to compute a spin configuration S and the corresponding energy Ti 

4. return S and H 
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FIG. 1: Ga,a grid graph. Dashed lines indicate negative edge 
weights (Color online). 
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remainder let F be the set of frustrated plaquettes in 
Gk,l and P the set of all plaquettes in Gk,l- 

Biechc et al. identify the frustrated plaquettes as 
vertices of a graph, Gf = (F,Ep) with F = {/ | 
/ is a frustrated plaquette in G} and Ep = F x F. Each 
edge e = (fi,fj) G Ep is assigned a weight w(e) equal 
to the sum of the absolute weights of the edges in Gk,l 
crossed by a minimum path connecting /; with fj. Fig- 
ure 3 shows the graph Gf for the grid graph of Figure 1. 
The underlying dual graph is shown in Figure 2. 

It is easy to see that minimizing the sum of the 
weights of unsatisfied edges connecting frustrated pla- 
quettes yields a spin configuration of minimum energy. 
The latter is achieved by determining a minimum- weight 
perfect matching in Gf- Finding a ground state is thus 
reduced to finding a minimum-weight perfect matching 
M of the graph Gf, and its energy is given as: 



H — JijSiSj 



<ho> 



= - £l J * 



- ]T \J ij \+2w(M) 



<»,i> 



For a detailed description of this method we refer to [1]. 



2. Limits of Bieche 's Approach 

The approach of Biechc ct al. is simple and intuitive, 
but comprised two major practical obstacles. First of all, 
in order to obtain the dual edge weights, one has to cal- 
culate shortest paths in Gd between all pairs of vertices. 



FIG. 2: Geometric dual graph Gd of the grid graph 
Gi,4 shown in Figure 1 which is seen translucent. Dark 
gray vertices represent frustrated plaquettes (assuming 
the trivial configuration <S° = {+1 | Vi G V}), and light 
gray vertices are unfrustrated plaquettes (Color online). 




FIG. 3: Graph Gf of the grid graph shown in Fig. 1. 
Continuous edges indicate distance 2 between vertices, 
and dotted edges indicate distance 1 (Color online). 



Although this can be done in time and space bounded 
by a polynomial in the size of the input, the calculations 
can take long in practice or require a large amount of 
memory. Equipped with the weights, one can construct 
the complete graph Gf of frustrated plaquettes. Treat- 
able system sizes are practically limited by the number 
of edges present in Gf- Assuming 32 bits for represent- 
ing an edge, one needs nearly 4 GB of memory just for 
representing Ep for a 300 x 300 grid graph, assuming 
50% of the plaquettes are frustrated. For a 400 x 400 
grid graph, almost 12 GB of memory are necessary, which 
goes beyond the hardware resources available in ordinary 
modern computers. 



4 



3. Simple Improvement for ± J -distributed Samples 

For ± J distributed instances, one can obtain the length 
of the shortest paths directly without shortest paths 
calculations. For this, we project the geometric dual 
graph Gd on the plane so that each vertex v is provided 
with definite coordinates (x v , y v ) preserving the distance 
function using rectilinear edges. Different vertices are 
assigned to different coordinates. Then the length of 
paths between i and j, ir — (xi, yi), a\, o<i . . . , a r , {xj,yj) 
(with <Ji = (xi,yi) = I G V) traversing only through 
the grid graph Gk.l without crossing its border, is 
given as the Manhattan distance c = \xi — Xj \ + \yt — 
yj\. This value is to be compared with the length 
of the path passing through the outer face o, tt — 
(xi,yi),a\, . . . ,ak,o,<Jk+i,- ■ -cr r , {xj,yj). The weight of 
this path is given as c Q = minja;^ y i: K — Xi, L — yi} + 
minjxj ,yj,K — Xj , L — yj}. The shortest path from i to j 
is the shorter of the two. In a half-torus graph, analogous 
calculations can be performed. 



4- The Approach of Barahona 

Barahona [8, 9] constructs the geometric dual graph 
that contains a vertex for each plaquette and edges 
in case the corresponding two plaquettes share an 
edge. Here the outer face is also interpreted as a 
plaquette. In formulas, Gd = (P,Ed) of Gk,l, 
where P = {p | p is a plaquette in Gk,l} U {o 
o is the outer face plaquette} and Ed = {e = (pi,Pj) | 
ypi,Pj € P, Pi H p j ^ 0}. Each dual edge is assigned 
a weight according to the absolute weight of the edge 
in Gk,l crossed by the dual edge. Vertices pi <G P are 
called odd if they represent a frustrated plaquette, oth- 
erwise even. 

Subsequently, the graph Gd is transformed into a 
graph G* . In order to do this, first every vertex pi G P 
with deg(pi) > 3 is expanded to (deg(pi) — 2) copies of 
degree 3. Any even vertex remains even, expanding an 
odd vertex makes one of its copies (arbitrarily) odd and 
the others even. From now on, one works with vertices 
of degree 3 only. Next, each vertex is transformed to 
a K3 subgraph: Each edge incident to an even vertex 
is replaced by an intermediate vertex and two edges. At 
most two new vertices are inserted for each edge connect- 
ing two even vertices. Original edges keep their weight, 
new edges obtain weight zero. For the details, we refer 
to [8, 9]. 

On G* a MWPM is computed. Any even vertex has 
an even number (including zero) of "outgoing" match- 
ing edges, however, any odd vertex has an odd count 
of those edges. After the matching is calculated, the 
afore expanded vertices are shrunken, and the remaining 
matching edges raise shortest paths connecting frustrated 
plaquettes. As the total length of the induced paths is 
minimal among all possible paths, the induced set of un- 
satisfied edges has minimum weight. Following Bieche, 



this corresponds to a configuration of minimum weight. 

5. Evaluation of Barahona 's Approach 

Barahona's transformation consists of two steps and is 
a bit more involved than the method of Bieche et al. For 
a quadratic L x L grid graph with free boundary condi- 
tions, G* has \V*\ « 12(£-l)[(L-l)+2]-12-6 0(W many 
vertices, where b dd — 3 • |{odd vertices} | is the number 
of odd vertices in the graph Gd- The graph G* is sparse 
as each vertex in G* has degree 3, \E*\ = §|V*|. Assum- 
ing 50% frustrated plaquettes, the number of vertices in- 
creases approximately by a factor of 10. Given that for 
bigger lattices this transformation needs less space than 
the one by Bieche et al., it is preferable to the former. 
However, in the next section we describe a method that 
works with an even smaller graph. 

6. The New Approach — Following in Kasteleyn's Footsteps 

In this section, we follow an idea first described by 
Kasteleyn [10, 20] and Fisher [21]. In these works, the 
goal was to calculate the configurational partition func- 
tion for dimer coverings on a lattice. The authors ex- 
ploited that the calculation of the partition function of 
the Ising model can be reduced to the number of ways in 
which a given number of edges can be selected to form 
closed polygons [22], i.e., a polygon configuration such 
that each lattice vertex has even degree of selected in- 
cident polygon edges. The latter can be computed as 
follows. First, one constructs a so-called cluster lattice 
graph which is generated by replacing each vertex of the 
lattice by a Kasteleyn city (a K4 subgraph) . Now the ex- 
panded graph is oriented such that the associated skew- 
symmetric Matrix D shows the property that |Pf(£))|, 
where Pf(D) denotes the Pfaffian of the skew-symmetric 
matrix, gives the number of dimer coverings. As there is a 
one-to-one correspondence between the number of poly- 
gon configurations and the number of dimer coverings, 
this method yields the generating function for polygon 
configurations and therefore the generating function for 
the Ising model. 

Thomas and Middlcton used Kasteleyn cities for cal- 
culating extended ground states on the torus in order 
to gain relevant information about the physics of spin 
glasses on toroidal lattices. Furthermore, they point out 
that the method yields an exact ground-state algorithm 
on planar lattices. 

A closely related approach was used later by Galluccio 
et al. to design an exact algorithm for the computation 
of the partition function for the Ising problem that runs 
in polynomial time for several models of interest [23- 
25], e.g., for two-dimensional toroidal lattices with ±J 
distribution. 

Here, we focus on planar grid graphs. The distribu- 
tion of the edge weights is arbitrary. Gk,l is transformed 
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FIG. 4: Expanded dual graph Gk 4 for the grid graph (^4,4 
(Color online). 



into a pseudo-dual graph Gk 4 as shown in Figure 4. For 
simplicity, we confine ourselves to quadratic grid graphs 
in the following, however all results can be easily trans- 
formed to general grid graphs. 

Formally, first the geometric dual of the grid graph 
Gl,l is constructed, then the outer face vertex is ex- 
panded into 2L + 1 (L in the half-torus case) copies, such 
that the resulting graph is an intermediate grid graph 
Gl+i.l+i (Gl,l+i)- The edge weights of Gl+i,l+i are 
set to the weights of the edges of Gl,l that are crossed by 
the edges of Gl+i,l+i- Edges that do not cross any other 
edge obtain weight zero. Next, each vertex of Gl+i,l+i 
is expanded to a K4 subgraph. Again, newly constructed 
edges receive weight zero. 

The transformation for the half-torus graph is done 
similarly, but the intermediate graph is a grid-half-torus 
graph Gl,l+i which is a grid with L — 1 additional edges. 
Edge weights are set as just described. Finally, all ver- 
tices are expanded to a K4, subgraph as before. We de- 
note by G 'inter the intermediate graph either for the un- 
derlying grid or half-torus graph. 

On the transformed graph Gk 4 we calculate a 
minimum- weight perfect matching M. 

The next step is to shrink all the ^-subgraphs back, 
resulting in the graph G inter- Also all copies of the 
outer face vertex are shrunken. Dealing again with the 
geometric dual graph of Gl,l, we take the subgraph 
Gs = {Q,5(Q)) of the geometric dual graph that con- 
sists only of dual edges that were matched, and all dual 
vertices with degree greater zero restricted to matched 
edges. This subgraph Gs is an Eulerian graph as each 
dual vertex is incident to an even number of matching 
edges. It is well known that there exists a one-to-one 
correspondence between Eulerian subgraphs in the dual 
graph and cuts in the original graph. So, Q defines a cut 
5(Q) in the graph Gl,l, cf. Fig. 5. 

In order to show the correctness of the transformation, 





FIG. 5: Backward transformation. A matching (thick lines) 
in Gk 4 induces a Eulerian subgraph in Gd and therefore a cut 
in the grid or half-torus graph. The vertex partition consists 
of the vertices with spin value +1 (+) in one partition and 
those with spin value —1 ( — ) in the other (Color online). 



we exploit that each ground state corresponds exactly to 
a MIN-CUT S(Q) in the grid graph Gl,l- A cut separates 
the vertex set into two disjoint sets W and V \ W . Ver- 
tices in the same partition get assigned the same spin 
value. Cut edges are those connecting a pair of vertices 
with different spin values. The Hamiltonian can be stated 



H — ^ JijSiSj 

<i,j> 



(2) 



<ho> 

We show that the edge set determined with the method 
described above corresponds to a minimum cut. 

w(M) = Y 

eeEDM 
w(e)=w(e) e£<5(Q) 

w{8{Q)) 
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As w(M) is the weight of a minimum- weight perfect 
matching, the weight of the subgraph G$ is minimum, 
and thus also the weight of the cut S(Q). 

As Kasteleyn's original method yields the complete 
partition function, one might ask why the algorithm was 
modified so that only the ground state is determined. 
The reason for this is twofold. First of all, minimum- 
weight perfect matchings can be computed in graphs with 
millions of vertices, provided they are sparse, which al- 
lows us to go to very large system sizes. Furthermore, 
the partition function does not encode the ground-state 
configuration itself but only its energy. 

7. Advantages of the New Approach 

The method is intuitive and its implementation is 
straightforward. We present some computational re- 
sults in Section V. For a quadratic L x x L x grid graph 
Glx,Lx we construct a graph G^ Ll with \V^ Ll \ = 

4 (L i + 1 ) ( L i + 1 ) vertices and 2 1 ,L 1 1 - 4 (L 1 + 1 ) edges . 
For an L 2 x L 2 ring graph Gl 2 ^ 2 the construction yields 

a graph G^ 2 with \V^' L2 \ = AL 2 {L 2 + 1) vertices and 

2|V^ 2,i2 | — 2L 2 — 4 edges. In any case, the resulting 
graphs are very sparse. More specifically, let us com- 
pare for an L x L grid the sizes of the graphs in which 
a minimum-weight perfect matching is computed. As- 
suming that around 50% of the original plaquettes are 
frustrated, the graph based on Kasteleyn cities contains 
only about one third of the number of vertices and one 
fourth of the number of edges contained in the graph 
constructed with Barahona's method. 

As the running time of the matching algorithm scales 
with the number of vertices and edges in the graph, the 
Kasteleyn construction is preferable to both the Bicchc 
et al. and the Barahona construction. 



Computing Domain Walls 

For the computation of domain walls, we follow the 
usual approach [26-28]. A ground state of the system 
is calculated, having energy Eq. To introduce a domain 
wall, the system is then usually perturbed by flipping all 
couplings along a row or a column in the lattice. The 
ground state for this new system is calculated, having 
energy Eq ert . The domain- wall energy for a given sample 
is then given by AE = \E^ ert — E \. We proceed as 
described in [29, 30] , by first determining a ground state 
of an EA spin glass with periodic boundary conditions 
in one direction, say along the x-axis. Then the signs 
of L edges in one column along the y-axis are flipped. 
The symmetric difference of these ground states yields a 
domain wall. 

With the help of linear programming [31, 32], one 
does not need to calculate the second ground state from 
scratch but can flip the weight of the specific L edges 



one by one, each followed by a reoptimization step. It is 
also possible to flip all signs at the same time and do a 
global reoptimization step. However, in practice for ±1 
distributed weights it often docs not pay off to do the 
reoptimization steps, and so we calculate both ground 
states from scratch. 



Modified Approach of Bieche et al. — Review of a 
Heuristic Method 

As argued above, the original approach of Bicche et 
al. suffers from an extensive memory usage. In order to 
overcome these limits, some modifications have been pro- 
posed that yield heuristic methods for low-energy states 
that are however not necessarily exact ground states. In 
these heuristics, a reduced graph G re d is used instead of 
a complete one. An approach often used is to introduce 
in G re d, only those edges with weight less than or equal 
to a fixed value c max (often c max = cJ max is chosen, with 
c = 4,5,6). For continuous spin systems, Weigel [33] re- 
cently suggested to introduce for each vertex only the k 
lightest edges. He used this cutoff-rule successfully for a 
matching routine embedded in a genetic algorithm. 

The reasoning behind this is that 'heavy'-weightcd 
edges are rarely contained in an optimum solution. This 
can assumed to be true for, e.g., ±1 distributed cou- 
plings and 50% negative weights, for which very often 
true ground states are reached in practice. 

Using the reduced graph G re d, Hartmann and Young 
determined high-quality heuristic ground states for L x L 
lattices with ±1 distribution. They could go up to L < 
480 [30]. Also using the heuristic variant based on the 
approach of Bieche et al. , Palmer and Adler report results 
for L < 1801 with the choice of c max = 6J max [34]. 

For different, especially smaller, percentage of negative 
weights, the quality of the heuristic decreases. This can 
be understood as follows. An edge (u,v) in the trans- 
formed graph is assigned a weight that equals the sum 
of the absolute weights of the edges in Gk,l crossed by 
a minimum path connecting u and v. u and v corre- 
spond to a pair of frustrated plaquettes. The latter can 
be assumed to be spread all over the system. In case 
of small p, the total number of frustrated plaquettes is 
small. Therefore, the weight of a minimum path connect- 
ing a pair of them can become large which can cause a 
heuristic with a limited value of c max to fail. 

Furthermore, for a different distribution of the cou- 
plings e.g., Gaussian couplings, this heuristic variant has 
to be used carefully, as good values for c max are not ev- 
ident. Certainly, removing heavy-weighted edges will re- 
sult in reduced graphs, but it is not clear beforehand 
which weights should be considered heavy. Applying dif- 
ferent cut-off rules, e.g. vertex-degree constraints, might 
be helpful as they were already used for thinning graphs, 
but suitable cut-off values depend still on the underlying 
distribution of the couplings. 

An experimental evaluation will be given in Section V. 
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In the next section, we describe how the correctness of 
heuristic ground states can be verified using linear pro- 
gramming. 



8. Checking Whether a Spin Configuration Actually 
Defines a Ground State 

Suppose we have a spin configuration at hand that has 
been computed by determining an optimum matching on 
the graph only containing 'light' edges and we want to 
test whether it is a correct ground state. 

First, we compute an optimum matching on the same 
reduced graph. In so-called pricing steps, we determine 
whether yet neglected edges exist that need to be taken 
into account in order to ensure correctness. In case such 
an edge is reported by pricing, it is introduced in the 
reduced graph. The process is iterated until no edge is 
returned any more. Pricing is a general feature in linear 
programming and combinatorial optimization. For more 
details, we refer to [31, 32, 35]. 

Pricing steps can be performed with Cook and Rohc's 
state-of-the-art Blossom IV code [36] by implementing 
small modifications. We give some computational results 
in Section V. 



V. COMPUTATIONAL RESULTS 

The method proposed in Section IV 6 can be used for 
any distribution of the couplings. Here, we focus on 
±1 distributed instances. The concentration p of anti- 
ferromagnetic bonds was set to 0.5. For the computa- 
tions we used Intel Xeon CPU with 2.3GHz and AMD 
Opteron Processor 248 with 2.2GHz, each with less than 
16GB RAM. The largest instances L > 1500 were done 
on Xeon processors with 16GB RAM. The physics anal- 
ysis done with our method will be presented in a forth- 
coming article [11]. 

Running times for different sizes of grid graphs (K = 
L = 100, 150, 250, 500, 700, 1000, 2000, and 3000) to- 
gether with the number of computed samples, are shown 
in Table I. For smaller grids with L < 500, we computed 
between 2*10 5 and 5*10 5 instances per size. For medium 
sized lattices 700 < L < 1500, we ran several thousand 
samples for each size, whereas for the largest systems of 
3000 2 spins we generated results for 157 samples. Lat- 
tices with L < 500 can be computed within less than 2 
minutes on average, whereas one ground-state determi- 
nation for the biggest size requires on average 24h CPU 
time on a single processor. Results of samples for spin 
glasses with periodic boundary conditions in one direc- 
tion are presented in Table II. The running times given 
in Tables I and II scale with L 3 - 50 ^ 5 ) for LxL spin glasses. 
In total we invested about 600 CPU days for our experi- 
ments. 

It turns out that free boundary samples are computa- 
tionally a bit more demanding than samples with peri- 







average time 


min-timc 


max-time 


nr. samples 


100 


0.67 


(±0.00) 


0.12 


3.47 


50925 


150 


2.03 


(±0.01) 


0.31 


16.80 


29750 


250 


9.70 


(±0.03) 


1.04 


121.18 


36800 


500 


109.62 


(±0.79) 


5.79 


1406.24 


20867 


700 


323.19 


(±4.54) 


18.96 


5233.04 


5499 


1000 


1200.33 


(±17.60) 


59.49 


9717.01 


3483 


1500 


5280.29 


(±111.79) 


288.58 


58036.33 


2330 


2000 


14524.34 


(±503.69) 


701.16 


117313.32 


942 


3000 


61166.70 


(±4920.19) 


3017.97 


316581.15 


157 





TABLE I: Running times and number of computed samples 
for different spin glass instance sizes with free boundary con- 
ditions and ±1 distributed couplings (p = 0.5). 





average time 


min-time 


max-time 


nr. samples 


400 


36.00 (±1.19) 


3.27 


712.67 


1400 


500 


88.31 (±2.91) 


5.65 


1239.81 


1900 


700 


319.38 (±3.53) 


17.46 


12712.34 


15847 


1000 


1129.03 (±27.20) 


49.87 


18577.58 


3000 



TABLE II: Running times and number of computed samples 
for different spin glass instance sizes with periodic boundary 
condition in one direction and ±1 distributed couplings (p = 
0.5). 



odic boundaries in one direction because the intermediate 
graphs G inter are larger. For other concentrations of anti- 
ferromagnetic bonds the presented data (running times, 
memory usage, etc.) are comparable, as our method, un- 
like the method of Bieche et al., does not depend on the 
concentration of anti-ferromagnetic bonds but only on 
the grid-graph size KL. 

Table III reflects the average over the maximal memory 
usage needed in the ground state calculations. The mem- 
ory usage roughly scales with L for Ising spin glasses 
with free boundary conditions and with L l e for Ising spin 
glasses with periodic boundaries in one direction. Both 
from the CPU time and from the necessary memory we 
conclude that the method is very fast. It needs consid- 
erably less memory than the commonly used method of 
Bieche et al., which in its heuristic variant allows to treat 
only smaller system sizes than reported here. However, 
we also note that a good statistics for system sizes beyond 
3000 2 would currently be hard to reach. 



Heuristic Ground States and Their Correction 

In this section, we explore the quality of the heuristic 
ground-state calculation using the method of Bieche ct 
al. for two-dimensional ±J Ising spin glasses with free 
boundary conditions. Within this we use the verification 
technique explained in Section IV 8. 

We consider planar LxL lattices with L = 164 and ±1 
distributed couplings. First we study these lattices with 
a concentration p = 0.5 of anti-ferromagnetic bonds. It 
turns out that out of 9912 computed samples 9 (0.091%) 
were wrong when using c max = 4. Thomas and Mid- 
dleton [12] stated 1.5% inexact solutions on toric sam- 
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1 fi^ 4 MR 


400 


24-5 6 MR 


500 


332.5 MB 


500 


321.6 MB 


700 


572.7 MB 


700 


544.6 MB 


1000 


994.6 MB 


1000 


993.7 MB 


1500 


2.052 GB 






2000 


3.568 GB 






3000 


7.832 GB 







TABLE III: Memory usage for different (±1, p = 0.5) sample 
sizes. 



pies with L < 128 and c = 8J max . We conclude that in 
practice the heuristic almost always returns true ground 
states if c max and p are suitably chosen. 

The overall average running time was 45.26 sec, com- 
paring an average of 82.97 sec. when pricing was nec- 
essary with 45.23 sec. without pricing. In our tests, 
one pricing step was always sufficient to correct a wrong 
ground state. Using the Kasteleyn approach, a ground 
state computation takes on average only around 1 second 
for this lattice size (on Xcon processors), as can be seen 
in Table V. 

In order to assess the influence of the cut-off parameter 
Cmax on the number of wrong results and the time to 
correct them, we varied p and c max for the ±1 164 x 164 
lattices (using the AMD Opteron processors). In Table 
IV we show results, always averaged over 100 instances. 

Several conclusions can be drawn from this experiment. 
Firstly, small values of c max lead to many wrong results. 
E.g., for c max = 3 up to 100% of the results were wrong, 
and the solutions have to be handled with care. Then, 
the quality of the heuristic increases with increasing c max . 
For large enough value, the solutions are very often cor- 
rect and only need a verification step in order to prove 
their correctness. 

Apart from this trend, the results suggest that the 
quality of the results highly depends on the chosen com- 
bination of c max and p. Clearly, for small percentage p 
one has to choose a higher cut-off value c max in order to 
generate high-quality solutions. E.g. forc max = 3, all 100 
results were wrong when the percentage of negative edges 
was chosen as p = 0.1, whereas 39 results were wrong for 
p = 0.5. This can also be seen by looking at different 
Cmax values but fixed p values. For c max < 6 and p = 0.1 
always at least one percent is wrong. This means that 
especially for smaller values of p, one has to make sure 
that the cut-off value is chosen big enough. The reason 
for this behavior is that the weight of a minimum path 
connecting a pair of frustrated plaquettes can become 
large for small p. Therefore, the minimum-weighted per- 
fect matching might contain heavy edges. c max has to be 
chosen large enough in order to ensure that these heavy 
edges are contained in the reduced graph. As argued 
before, the necessity of having to choose high cut-off val- 



% edges w. 


average 


time [sccj 


pricing 


avg. nr. of 


nr. of 


w(e) = -1 






time [sec] 


pricing steps 


wrong 
results 


Cmax = 3 


10 


16.91 


(±0.70) 


9.10 


2.03 


100 


20 


51.59 


(±1.13) 


18.30 


1.11 


93 


30 


56.87 


(±1.57) 


22.42 


1.00 


49 


40 


56.95 


(±1.59) 


24.01 


1.00 


36 


50 


58.72 


(±1.69) 


25.10 


1.00 


39 






t 


max = 4 






10 


15.60 


(±0.35) 


4.76 


1.09 


77 


20 


37.37 


(±0.74) 


17.99 


1.00 


6 


30 


49.93 


(±0.67) 


0.00 


0.00 





40 


53.74 


(±0.72) 


0.00 


0.00 





50 


54.02 


(±0.76) 


0.00 


0.00 









c 


max = 5 






10 


13.04 


(±0.30) 


4.64 


1.00 


18 


20 


37.96 


(±0.55) 


0.00 


0.00 





30 


52.41 


(±0.93) 


0.00 


0.00 





40 


56.79 


(±1.00) 


0.00 


0.00 





50 


59.30 


(±1.01) 


0.00 


0.00 











max — 6 






10 


13.11 


(±0.21) 


5.04 


1.00 


1 


20 


42.52 


(±0.71) 


0.00 


0.00 





30 


59.05 


(±1.30) 


0.00 


0.00 





40 


65.14 


(±1.44) 


0.00 


0.00 





50 


62.26 


(±1.16) 


0.00 


0.00 









c 


max = 7 






10 


14.06 


(±0.23) 


0.00 


0.00 





20 


46.70 


(±0.89) 


0.00 


0.00 





30 


66.18 


(±1.65) 


0.00 


0.00 





40 


73.18 


(±1.95) 


0.00 


0.00 





50 


69.69 


(±1.52) 


0.00 


0.00 











max — 8 






10 


15.01 


(±0.29) 


0.00 


0.00 





20 


53.10 


(±1.32) 


0.00 


0.00 





30 


74.54 


(±2.15) 


0.00 


0.00 





40 


83.20 


(±2.32) 


0.00 


0.00 





50 


82.06 


(±2.13) 


0.00 


0.00 






TABLE IV: Different c max values used within the heuristic 
version of the Bieche et al. approach for each 100 ±1 164 x 164 
grid graphs with various concentration of anti-ferromagnetic 
bonds. 



ues can lead to memory problems as the edge density 
of the generated reduced graphs increases. These diffi- 
culties can be avoided by using the approach based on 
Kasteleyn cities. 

From Table IV we see that pricing takes only negligible 
running time. As a conclusion, if the Bieche et al. algo- 
rithm is used to determine ground-state properties, it is 
advantageous to do the pricing steps, too, independent 
of the size of the reduced graph G re d- However, despite 
the fact that pricing is very fast, the heuristic together 
with the verification is still considerably slower than the 
method proposed here based on Kasteleyn cities, cf. Ta- 
ble V. 

It is also interesting to assess the quality of the heuris- 
tic for different distribution of the couplings, e.g. for 
Gaussian distributed couplings. We study grid graphs 
Gl,l with L = 50, 100, 150. The grid graph size L was 
limited due to the fact that the generation of the graphs 
Gf takes a long time. This is explainable as one is in 
the need of all-pair shortest path calculations for the set 
of frustrated vertices on the dual graph. The latter is 
a polynomially solvable problem, however the computa- 
tions can take long. In our tests, computing the short- 
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10 


20 


30 


40 


50 


164^ 


0.92 (±0.005) 


1.28 (±0.005) 


1.32 (±0.007) 


1.17 (±0.007) 


1.11 (±0.009) 



TABLE V: Average running times (in sec.) over 10.000 instances, using the method based on Kasteleyn cities on ±1 164 x 164 
grids with various concentration of anti-ferromagnetic bonds. 



est paths usually takes much longer than computing the 
ground states itself. In Table VI the running times for 
the matching on the different reduced graphs can be seen. 

As it is not clear beforehand how big the cut-off pa- 
rameter c max should be, we let ourselves be guided by the 
corresponding values used in instances with bimodal dis- 
tribution. More specifically, for some value of c max (we 
used c max = 5) we compute the average percentage of 
'light' edges from Gf that are contained in the reduced 
graph G rec i in instances with bimodal distribution. For 
L = 50, we find that in ±1 distributed instances with 
p = 0.5 on average the 8% lightest edges are used (this 
percentage reduces to an average of 2% for L — 100 and 
to 1% for L = 150). 

For Gaussian distributed instances, we build the re- 
duced graphs with the same percentages of light edges, 
i.e., 8% (2% or 1%) light edges for L = 50 (L = 100, 
L = 150 resp.) Results are shown in Table VI. First 
of all, it turns out that for small grid graphs many re- 
sults are wrong. For L = 50, about 12.5% of the calcu- 
lated instances do not find correct ground states, and a 
higher cut-off value has to be used. Considering larger 
grid graphs the situation looks similar. 47% (52%) in- 
stances computed the wrong ground state for L = 100 
(L = 150). 

For L — 50, we have to go up to a percentage of 16% 
lightest edges (corresponding to c max ~ 8), for L = 100 
to 4% (corresponding to c max s=s 7) and for L = 150 
we have to take into account the 2.5% lightest edges, 
corresponding to c max w 8, in order to ensure correctness 
of the results for our test data. Again, when comparing 
the running times for 150 2 lattices, it is by roughly a 
factor of 25 faster to use the Kasteleyn city approach 
instead of a heuristic Bieche method. 

As the performance of the matching routine scales with 
the graph sizes, more specifically with the number of ver- 
tices and edges, we study the reduced graphs with respect 
to these two entities. Usually, the graphs for computing 
the matchings are dense. This is especially true for small 
grid graphs (L — 50), where 16% of the light edges are 
needed. Increasing the grid graphs leads to a consider- 
able decrease of needed light edges, and decreasing den- 
sity. (The density of a graph with n vertices is defined as 
the number of its edges divided by the number of edges 
of the complete graph K n .) Nevertheless, in our exper- 
iments the grid graphs with L = 100 (L = 150) contain 
on average \V\ = 4901 ± 5 (|V| = 11103 ± 8) vertices. 
Taking 4% (2.5%) of all edges means in absolute num- 
bers 4.803(9) * 10 5 (1.541(2) * 10 6 , resp.) edges, which are 
large and dense graphs. This has to be compared with 



the graphs used within the Kasteleyn approach. Here we 
have for L = 100 (L = 150) graphs with \V\ = 40804 
(\V\ = 91204) and 8.120 * 10 4 (1.818 * 10 5 ) edges. In 
fact, the graphs we deal with have more vertices but 
are considerably sparser. The Blossom IV routine can 
compute matchings in very large graphs, provided their 
density is low. This fact is also reflected in the better 
running times for the Kasteleyn approach presented in 
this section. Comparing with the situation for instances 
with bimodal distributed edge weights, we see from Ta- 
ble IV that a value c max = 4 in the case of p = 0.5 
yields good results. For L = 100 (L = 150) we have 
|Vj = 4901±4 (|Vj = 11114±7) with \E\ = 1.58(7)* 10 s 
(\E\ = 3.66(1) * 10 5 ), which arc about 1.3% (0.6%) edges 
of the complete graph. Thus, the graphs with bimodal 
distribution can be chosen sparser than in the Gaussian 
case (for c max = 4). However, these graphs are usually 
denser than the graphs used in the Kasteleyn approach. 



VI. CONCLUSIONS 

We presented a simple algorithm (Section IV 6) based 
on Kasteleyn cities. The algorithmic foundations of this 
method date back to the work of Kasteleyn [10] from the 
1960s in which he computed the complete partition func- 
tion for the Ising model. Using this approach, we can 
compute exact ground states for two-dimensional planar 
Ising spin-glass instances. The method is easy to imple- 
ment, fast and has only limited memory requirements. 
According to our knowledge, the treatable system sizes 
are considerably bigger than the ones computed earlier 
and are always provably exact. Thomas and Middleton 
[12] used Kasteleyn cities for studying extended ground 
states. Furthermore, they state that the method can also 
be used for determining exact ground states of planar 
graphs. 

We evaluated different established exact methods and 
compared them with respect to running time and mem- 
ory requirements. It turned out that the approach pre- 
sented here is both considerably faster and needs less 
memory than the methods proposed earlier. We showed 
how heuristically computed ground states can be verified 
or corrected fast using mathematical optimization. How- 
ever, the method based on Kasteleyn cities still outper- 
forms this approach. Finally, we evaluated the solution 
quality of heuristic variants of the Bieche et al. approach. 

In the future, we will make our program available in 
public domain via the Cologne Spin Glass Server that 
can be found at 
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L 


Cmax (±1 case) 


% lightest edges (Gassian case) 


% wrong results 


matching time [sec] 


50 


5 


8.0 


12.5 


0.11 (± 0.01) 




6 


12.0 


2.0 


0.18 (± 0.01) 




7 


15.0 


2.0 


0.23 (± 0.01) 




8 


16.0 


0.0 


0.24 (± 0.01) 


100 


5 


2.0 


47.0 


2.68 (± 0.20) 




6 


3.0 


2.0 


7.13 (± 0.50) 




7 


4.0 


0.0 


11.27 (± 0.80) 




8 


5.0 


0.0 


23.20 (± 1.60) 


150 


5 


1.0 


52.0 


34.60 (± 2.21) 




6 


1.4 


7.0 


49.34 (± 2.97) 




7 


1.9 


0.0 


66.70 (± 4.00) 




8 


2.5 


0.0 


85.18 (± 5.14) 



TABLE VI: Size of reduced graph G rec i when using Gaussian distributed couplings compared to the achieved quality by the 
heuristic variant of the Bieche et al. approach. 



http : //cophy . informatik.uni-koeln.de/. 
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